﻿module Bayes



open System
open System.IO
open MathNet.Numerics
open MathNet.Numerics.FSharp
open MathNet.Numerics.LinearAlgebra.Double
open System.Diagnostics
open System.Collections.Generic

//struttura dati di appoggio
type GaussianPoint =
  struct 
    val Point   : DenseVector
    val Label   : float
    new (point, label) = { Point = point; Label = label }
  end

//funzioni ausiliarie
let fi (y:float) (points : List<GaussianPoint>) =
  points |> Seq.filter(fun p -> p.Label = y)  
         |> Seq.length 
         |> (fun sum -> sum / points.Count)


let mu (y:float) (points : List<GaussianPoint>) =
  let empty_vector = DenseVector(points.[0].Point.Count)
  let y_points = points |> Seq.filter(fun p -> p.Label = y) 

  y_points |> Seq.fold(fun k p -> k + p.Point) empty_vector
           |> (fun sum -> sum / (y_points |> Seq.length |> float))

let p_k (input : List<GaussianPoint>) (k : float) =
  let k_points = List(input |> Seq.filter(fun p -> p.Label = k))
  (float k_points.Count) / (float input.Count)


let sigma2_ki (k:float) (i:int) (points : List<GaussianPoint>) (mu_list:List<DenseVector>) =
  let k_points = List(points |> Seq.filter(fun p -> p.Label = k))
  k_points  |> Seq.mapi(fun i p -> (p,i))
            |> Seq.fold(fun res (p,j) -> let vector = p.Point.[i] - mu_list.[int k].[i]
                                         let squared_vector = vector * vector
                                         squared_vector + res) 0.
            |> (fun sum -> sum / (float k_points.Count))

let prob_bayes_ki (observed_point : DenseVector) (k : float) (i : int) (points : List<GaussianPoint>) (mu_list:List<DenseVector>) =
  let sigma2 = sigma2_ki k i points mu_list
  (1. / Math.Sqrt(2.0 * Math.PI * sigma2)) * Math.Exp(- (Math.Pow(observed_point.[i] - mu_list.[int k].[i], 2.)) /(2.0 * sigma2) )
  
 
let compute_bayes (training_set : List<GaussianPoint>) (y_list : seq<float>) =
  let y_list = List(y_list |> Seq.map(fun k -> k, p_k training_set k))
  let mu_list = List(seq{ for y in y_list do yield mu (fst y) training_set})
  seq{
    for i in [0..149] do
    
      let observed_point : DenseVector = training_set.[i].Point
      let p_xy =
        List(seq{for (k, p_k) in y_list do
                  yield 
                    List(seq{for i in [0..(observed_point.Count-1)] do yield (prob_bayes_ki observed_point k i training_set mu_list) * p_k }) |> Seq.fold(fun res p -> p * res) 1.
        })

      //stampa del risultato dell'osservazione
      let computer_label = 
          if p_xy.[0] > p_xy.[1] && p_xy.[0] > p_xy.[2] then
            Console.WriteLine("0")
            0.
          elif p_xy.[1] > p_xy.[0] && p_xy.[1] > p_xy.[2] then
            Console.WriteLine("1")
            1.
          else
            Console.WriteLine("2")
            2.

      yield (training_set.[i].Point, training_set.[i].Label, computer_label)
  }
